pacman::p_load(tidyverse,data.table,broom,ggpubr,latex2exp,ggpubr,sf,raster)
rm(list=ls())
sf_use_s2(FALSE)
################################################################################


########################## Time series

df = read_csv(paste0('../../data/landuse/clean/cropland_southamerica_nation.csv.gz'),show_col_types = FALSE)
df$area_planted_ha =  df$area_planted_ha/1000000
vartitle = 'A. Area planted in Argentina and Brazil (millions of ha)'
crop_list <- c("soybean","maize", "wheat","sugarcane","rice")
dfs <- df %>% filter(crop %in% crop_list)
dfs$crop <- factor(dfs$crop, levels = crop_list)
dfs <- dfs[,c('year','crop','area_planted_ha')]
setnames(dfs, old='area_planted_ha',new='variable_name')
setDT(dfs)

figure2a = ggplot(data = dfs[year<=2018,],  aes(x=year, y=variable_name, color=crop) ) +
  geom_line(linewidth=0.25) +
  geom_point(aes(shape=crop), size=2) +
  labs(title=vartitle) +
  theme_minimal() +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size=8, color='black'),
        axis.text.y = element_text(size=8, color='black'),
        plot.title = element_text(size=9),
        legend.position = "inside",
        legend.position.inside = c(0.11,0.7),
        legend.title = element_blank(),
        legend.text = element_text(size = 8),
        legend.key.width = unit(0.75,"cm"),
        legend.key.height = unit(0.5,"cm"),
        aspect.ratio = 0.35,
        panel.grid.minor = element_blank()
  ) +
  scale_color_manual(values=c("#D55E00", "#009E73", "#E69F00","#CC79A7","#0072B2"))+
  scale_x_continuous(breaks=seq(1990, 2018, 4))

figure2a_bw <- figure2a + scale_color_grey(start=0, end=0.5)


########################## Cross section

for(country in c('brazil','argentina')){
  
  df_l = fread(paste0('../../data/landuse/clean/landuse_',country,'_county.csv.gz'))
  df_u = fread(paste0('../../data/landuse/clean/geographicunits_',country,'.csv.gz'))
  
  if (country=='brazil'){
    start_year = '1995'
    end_year = '2017'
    df_u = df_u[, list(county_id, amc_id, microregion_id, mesoregion_id, state_id, region, land_area_km2)]
  } else {
    start_year = '2002'
    end_year = '2018'
    df_u = df_u[, list(county_id, state_id, region, land_area_km2)]
  }
  df = merge(df_l, df_u, by='county_id')
  
  df = df[, area_id := county_id]
  if (country=='brazil'){
    df = df[, c('county_id','amc_id','microregion_id','mesoregion_id', 'state_id') := list(NULL,NULL,NULL,NULL,NULL)]
  } else {
    df = df[, c('county_id','state_id') := list(NULL,NULL)]
  }
  df = df[, lapply(.SD, sum, na.rm=TRUE), by=list(year, region, area_id)]
  
  if (country=='brazil'){
    df = df[, c('farmland','pasture','forest') := list(
      crops_permanent + crops_seasonal2 + pasture_natural + pasture_planted + forest_natural + forest_planted + other,
      pasture_natural+pasture_planted, forest_natural+forest_planted)]
  } else {
    df = df[, c('farmland','forest') := list(
      crops_permanent + crops_seasonal2 + crops_unaccounted + pasture + forage_seasonal + forage_permanent + forest_natural + forest_planted + not_used + not_usable + trails_and_parks + unaccounted,
      forest_natural+forest_planted)]
  }
  df = df[, c('soybean_share','pasture_share','forest_share'):=list(soybean/(farmland), pasture/(farmland), forest/(farmland))]
  
  df_x = df[year==start_year]
  df_y = df[year==end_year]
  dff = merge(df_x, df_y, by=c('region', 'area_id'))
  dff = dff[,c('delta_soybean_share','delta_pasture_share','delta_forest_share') := list(soybean_share.y-soybean_share.x, pasture_share.y-pasture_share.x, forest_share.y-forest_share.x)]
  dff[soybean.x==0 & soybean.y==0, delta_soybean_pc := 0]
  dff[pasture.x==0 & pasture.x==0, delta_pasture_pc := 0]
  dff[forest.x==0 & forest.y==0, delta_forest_pc := 0]
  dff = dff[, list(region, area_id, delta_soybean_share, delta_pasture_share, delta_forest_share)]
  
  if(country=='brazil'){dff_br = dff}else{dff_ar = dff}
  
}

#Maps
shpfile = sf::read_sf("../../data/shapefiles/clean/argentinabrazil_county.shp")
shpfile <- st_crop(shpfile, extent(-85,-33, -65, 11.75))
shpfile <- st_simplify(shpfile, dTolerance = 0.1, preserveTopology = T)
shpfile_ctry = sf::read_sf("../../data/shapefiles/clean/restofsouthamerica_country.shp")

for (fill_var_name in c('delta_soybean_share','delta_pasture_share','delta_forest_share') ){
  
  df = setDT(rbind(dff_ar, dff_br))
  setnames(df, fill_var_name, "fill_var")
  df = df %>% drop_na(fill_var)
  df = df[fill_var!=0,]
  df = mutate(df, fill_var = ifelse(fill_var > 0.1, 0.1, fill_var))
  df = mutate(df, fill_var = ifelse(fill_var < -0.1, -0.1, fill_var))
  dfm = merge(shpfile, df, by.x="cty_id", by.y="area_id")
  limit <- max(abs(dfm$fill_var))*c(-1, 1)
  
  if(fill_var_name=='delta_soybean_share'){fill_var_lab = TeX("B. $\\Delta$ soybean share of land")}
  if(fill_var_name=='delta_pasture_share'){fill_var_lab = TeX("C. $\\Delta$ pasture share of land")}
  if(fill_var_name=='delta_forest_share'){fill_var_lab = TeX("D. $\\Delta$ forest share of land")}
  
  ggp = ggplot() +
    geom_sf(data=shpfile_ctry, fill=NA, color = 'black', linewidth = 0.05) +
    geom_sf(data=shpfile, fill=NA, color = 'gray', linewidth = 0.05) +
    geom_sf(data=dfm, aes(fill=as.vector(fill_var), color = after_scale(fill)), linewidth = 0.05) +
    scale_fill_fermenter(palette = "RdYlGn", direction = 1, limit=limit) +
    theme_void() +
    labs(title = fill_var_lab, fill = '') + 
    theme(legend.position = "bottom",  legend.direction = "horizontal", 
          legend.text = element_text(size= 7),
          legend.key.width = unit(0.6,"cm"),
          legend.key.height = unit(0.2,"cm"),
          plot.title = element_text(size = 9)) 
  
  ggp_bw = ggplot() +
    geom_sf(data=shpfile_ctry, fill=NA, color = 'black', linewidth = 0.05) +
    geom_sf(data=shpfile, fill=NA, color = 'gray', linewidth = 0.05) +
    geom_sf(data=dfm, aes(fill=as.vector(fill_var), color = after_scale(fill)), linewidth = 0.05) +
    scale_fill_fermenter(palette = "Greys", direction = 1, limit=limit) +
    theme_void() +
    labs(title = fill_var_lab, fill = '') + 
    theme(legend.position = "bottom",  legend.direction = "horizontal", 
          legend.text = element_text(size= 7),
          legend.key.width = unit(0.6,"cm"),
          legend.key.height = unit(0.2,"cm"),
          plot.title = element_text(size = 9)) 
  
  if(fill_var_name=='delta_soybean_share'){ggp_s=ggp;ggp_s_bw=ggp_bw}
  if(fill_var_name=='delta_pasture_share'){ggp_p=ggp;ggp_p_bw=ggp_bw}
  if(fill_var_name=='delta_forest_share'){ggp_f=ggp;ggp_f_bw=ggp_bw}
}

figure2b = ggarrange(ggp_s,ggp_p,ggp_f, nrow=1,ncol=3)
figure2b_bw = ggarrange(ggp_s_bw,ggp_p_bw,ggp_f_bw, nrow=1,ncol=3)

figure2 = ggarrange(figure2a,figure2b, nrow=2,ncol=1)
ggsave(paste0("../../output/figures/figure2.pdf"), width = 5.5, height = 5, units='in')
ggsave(paste0("../../output/figures/figure2.eps"), width = 5.5, height = 5, units='in')

figure2_bw = ggarrange(figure2a_bw,figure2b_bw, nrow=2,ncol=1)
ggsave(paste0("../../output/figures/figure2_bw.pdf"), width = 5.5, height = 5, units='in')
ggsave(paste0("../../output/figures/figure2_bw.eps"), width = 5.5, height = 5, units='in')



########################## Exports

#World imports
df_main = fread(paste0('../../data/trade/clean/imports_from_world_iv.csv.gz'))
df = df_main[item=='SOYBEANS',]
df = df[unit=='TONNES',]
df = df[destination!="CHINA",]
df$imports_from_argentina_brazil = df$imports_from_argentina + df$imports_from_brazil
df = df[,c('year','imports_from_world','imports_from_argentina','imports_from_brazil','imports_from_argentina_brazil')][, lapply(.SD, sum, na.rm=TRUE), by=list(year)]
setnames(df, new=c('year','World','World (from Argentina)','World (from Brazil)','World (from Argentina and Brazil)'))
df_w = melt(df, measure.vars = c('World','World (from Argentina)','World (from Brazil)','World (from Argentina and Brazil)'),
            variable.name = "destination", value.name = "tonnes_imported")
df_w$destination_group = 'World'

#China imports
df = df_main[item=='SOYBEANS',]
df = df[unit=='TONNES',]
df = df[destination=="CHINA (MAINLAND)",]
df$imports_from_argentina_brazil = df$imports_from_argentina + df$imports_from_brazil
df = df[,c('year','imports_from_world','imports_from_argentina','imports_from_brazil','imports_from_argentina_brazil')]
setnames(df, new=c('year','China','China (from Argentina)','China (from Brazil)','China (from Argentina and Brazil)'))
df_c = melt(df, measure.vars = c('China','China (from Argentina)','China (from Brazil)','China (from Argentina and Brazil)'),
            variable.name = "destination", value.name = "tonnes_imported")
df_c$destination_group = 'China'

df = rbind(df_c,df_w)
ggplot(data = df[destination %in% c('World','China','China (from Argentina and Brazil)'),], aes(x=year, y=tonnes_imported/1000000, color=destination, linetype=destination, shape=destination)) +
  geom_line(linewidth=0.25) +
  geom_point(aes(color=destination, shape=destination), size=3) +
  labs(title=paste0('Soybean imports (million metric tonnes)')) +
  theme_minimal() +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.title = element_blank(),
        legend.position = "inside",
        legend.position.inside = c(0.25,0.8),
        title = element_text(size = 9),
        aspect.ratio = 0.5,
        panel.grid.minor = element_blank())+
  scale_x_continuous(breaks=seq(1990, 2018, 2))+
  scale_color_manual(values = c("darkred","darkred","black"))+
  scale_linetype_manual(values = c("solid","dashed","solid"))
ggsave(paste0("../../output/figures/appendix/figure9_chinademand.pdf"), width = 7, height = 4)
ggsave(paste0("../../output/figures/appendix/figure9_chinademand.eps"), width = 7, height = 4)

